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An impact of particles' roughness on the self-diffusion 
coefficient D in granular gases is investigated. For a sim- 
plified collision model where the normal, e, and tangen- 
tial, /3, restitution coefficients are assumed to be constant 
we develop an analytical theory for the diffusion coeffi- 
cient, which takes into account non-Maxwellain form of 
the velocity-angular velocity distribution function. We 
perform molecular dynamics simulations for a gas in a 
homogeneous cooling state and study the dependence of 
the self-diffusion coefficient on e and /3. Our theoretical 
results are in a good agreement with the simulation data. 



1 Introduction 

Among numerous contributions by Isaac Goldhirsch 
to the theory of granular fluids |10) are his pioneering 
works on the hydrodynamic of gas of particles with 
the rotational degrees of freedom [IT]. We wish to 
dedicate this article, addressed to granular gases of 
rough particles, to the memory of Isaac Goldhirsch. 

Granular fluids are systems composed of a large 
number of macroscopic particles which suffer dissipa- 
tive collisions. In many respects they behave like or- 
dinary molecular fluids and may be described by the 
standard tools of kinetic theory and hydrodynamics 
[10[ [5j [14] . Among numerous phenomena, common 
to molecular and granular fluids, are Brownian mo- 
tion, diffusion and self-diffusion H H E3 H H 01] ; 
in the latter case tracers are identical to surrounding 
particles. The self-diffusion coefficient D has a mi- 
croscopic and macroscopic (thermodynamic) mean- 



ing, e.g. [HIS]: Microscopically, it characterizes the 
dependence on time of the mean-square displacement 
of a grain (R 2 (t)}: 



(R 2 (t)) 



6 J D(t')dt' , 
o 



(1) 



while macroscopically, it relates the macroscopic flux 
of tracers J s (r, t ) (the index s refers to tracers) to 
its concentration gradient Vn s (r, i): 



J 8 {r,t) = ~D{t)Vn s (r,t) 



(2) 



With the continuity equation " n, g' t ' z ' + V J s (r,t ) =0 



dn s (r,t) 

the local tracers density n s (r, t ) obeys the diffusion 
equation: 



dn s (r ) 
dt 



DV 2 n s (r ) 



(3) 



If the state of a granular gas is not stationary, like, 
e.g. a homogeneous cooling state, the diffusion coef- 
ficient generally depends on time. 

In the previous studies the self-diffusion coefficient 
has been calculated only for smooth particles, that 
is, for particles which do not exchange upon collisions 
the energy of their rotational motion [51 l4l ITJ1 IfJl 121 fTTj] . 
The macroscopic nature of granular particles, how- 
ever, implies friction between their surfaces and hence 
the rotational degrees of freedom are unavoidably in- 
volved in grains dynamics. 

In the present paper we analyze the impact of par- 
ticles roughness on the self-diffusion in granular gases 
- rarified systems, where the volume of the solid frac- 
tion is much smaller than the total volume. We cal- 
culate the self-diffusion coefficient as a function of 
normal and tangential restitution coefficients, which 
we assume to be constant. We develop an analytical 
theory for this kinetic coefficient and perform molec- 
ular dynamics (MD) simulations. The rest of the ar- 
ticle is organized as follows. In Sect. [2] we formulate 
the model, introduce the necessary notations and give 
the detailed description of the theoretical approach. 
In Sect. [3] we present the results of MD simulations 
and compare the numerical data with the predictions 
of our theory. In the last Sect. [4] we summarize our 
findings. 



2 Calculation of the 
diffusion coefficient 

2.1 Collision rules 



The dissipative collisions of grains are characterized 
by two quantities, e and (3 - the normal and tangen- 
tial restitution coefficients. These coefficients relate 
respectively the normal and tangential components 
of the relative velocity between surfaces of colliding 
grains, g — V12 + § (e x (aJi +uj 2 )), before (unprimed 
quantities) and after (primed quantities) a collision 

mm- 



self- 2.2 Boltzmann equation 

To compute the self-diffusion coefficient in granular 
gases two main approaches may be exploited - Green- 
Kubo approach, based on the time correlation func- 
tions of a dynamical variable (a particle velocity for 
the case of D) and Chapman-Enskog method, where 
the evolution of velocity distribution function is ana- 
lyzed. It may be shown that these two are equivalent 
for the case of constant normal restitution coefficient 
[5J |5]; here we adopt the latter approach. First we 
consider the force-free granular gas in a homogeneous 
cooling state. We start from the Boltzmann equation 
for the distribution function of tracers f s (v,u,r,t), 



(g n r 
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Here V12 = V\—V2 is the relative translational velocity 
of particles, uj\ and L02 are their rotational velocities, 
the unit vector e is directed along the inter-center vec- 
tor f*i2 at the collision instant and a is the diameter 
of particles. The normal and tangential components 
of g read respectively, g n = (g ■ e) e and g = g — g n . 
For macroscopic particles the normal restitution coef- 
ficient e can vary from (completely inelastic impact) 
to 1 (completely elastic collision^ while the tangen- 
tial coefficient j3 ranges from —1 (absolutely smooth 
particles) to +1 (absolutely rough particles). From 
the conservation laws and Eqs. Q follow the after- 
collision velocities and angular velocities of particles 
in terms of the pre-collision ones, e.g. [51 |2"T1 IT2"] : 



«1,2 = v l,2 T 
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Here 77 



2(1+9) 



characterizes friction, q = AI/ma 2 



with I and m being respectively the moment of in- 
ertia and mass of particles. Although the normal 
and tangential restitution coefficients generally de- 
pend on the translational and rotational velocities 
of colliding particles and an impact geometry, e.g. 
[201 [H E] , here we assume for simplicity that 
the restitution coefficients arc constants. 



^^For nano-particles the normal restitution coefficient can 
attain negative values [17] . 
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where f(v, u, f, t) is the velocity distribution function 
of the gas particles, 52 (c) is the contact value of the 
pair correlation function, which accounts for the in- 
creasing collision frequency due to the excluded vol- 
ume of grains and I(f,f s ) abbreviates the collision 
integral El US]: 

I{f,fs) = o 2 J dv 2 J d,Lu 2 JdeO(-vi2 ■ e) |ui 2 • e \ x 



■U (yl', CSi') / (v 2 ", £ 2 ") - f s (v u dJi) / (v 2 , cD 2 ) 



e 2 /3 2 



(7) 



Here v{' 2 and u" 2 denote the pre-collision velocities 
and angular velocities of particles for the inverse colli- 
sion, i.e., for the collision which ends up with vi t 2 and 
UJ1.2, the factor (— vi 2 • e ) selects only approaching 
particles, and \vi2 ■ e \ gives the length of the colli- 
sion cylinder. We assume that the concentration of 
tracers 



n s (r,t)= / dvi / duif s (vi,ui,r,t) 



(8) 



is much smaller than the number density of the gas 
particles n(r,t) = n = const, so that tracers do not 
affect the velocity distribution function / for the gas, 
which is the solution of the Boltzmann equation 



8_ 
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f(v,u,t)=g 2 (a) /(/,/) 



(9) 
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for the homogeneous cooling state. The collision in- 
tegral in the last equation is given by Eq. ([7| with 
the distribution function of tracers f s substituted by 
that for the gas particles /. 

Solving the Boltzmann equation ^ for 
f s (v, uj, r, t), which is also known as Boltzmann- 
Lorentz equation, one can write the diffusion 
flux 

J s (r f ,t) = J dv J dujvf s (v,uj,f,t) , (10) 

and then, using the macroscopic equation find 
the diffusion coefficient as the coefficient at the con- 
centration gradient Vn s {r,t). 

In a homogeneous cooling state a density of the 
system remains uniform, while the translational T{t) 
and rotational R(t) granular temperatures, defined 
through the average energies of translational and ro- 
tational motion, respectively, 



—nT — f dvduj f, 

2 / 2 



: nR = f duduj / 



(11) 

decay with time as dT/dt = —(T and dR/dt = 
— CrT, where £ and £r are the translational and ro- 
tational cooling rates, e.g. ETJ [7J [18] . 

The distribution function / (v, uj, t) is close to the 
Maxwell distribution and the deviations from the 
Max-wellian are twofold. Firstly, the translational 
and angular velocity distributions slightly differ from 
the Gaussian distribution, which is accounted by the 
expansion of / in Sonine polynomials series with re- 
spect to v 2 and uj 2 . Secondly, a linear velocity and 
spin of particles are correlated; this is quantified us- 
ing an expansion in Legendre polynomials series with 
respect to 6 = vuj - the angle between the velocity 
v and angular velocity uj of a particle. Using the 
reduced quantities c = v/vt and w — lo/lot with 
vt = y/2T/m and u>t = y/2R/I, the distribution 
function reads 



f(v,uj,t) 



-K 2 v\ Uj\ 



~ w2 [l + a w s[% (c 2 ) + 

a m S\y 2 (w 2 ) + a 11 s[ 1 / 2 (c 2 )S { 1 1 / , 2 (w 2 ) + 

bc 2 w 2 P 2 (cos6)}, (12) 



The Sonine and Legendre polynomials in Eq. ( 12 1 are 

S^hx) = | - x, Sfhx) = |(15 - 20a; - Ax 2 ) and 

Pzix) = | (x 2 — I)- The distribution function at- 
tains the scaling form ( [l2| ) after some relaxation time. 
At this stage, which will be addressed below, the co- 
efficients O2o, a iij a 02 and b are constants. These co- 
efficients have been analyzed in detail in Refs. [T8"HT]. 
They are complicated nonlinear functions of e, f3 and 
r = R/T, which are too cumbersome to be given here. 

The temperature ratio r also attains a steady state 
value in the scaling regime and the cooling rates be- 
come equal, ( = (r [HI E] • The steady-state value of 
r = R/T may be obtained as an appropriate root of 
the eighth-order algebraic equation [TH]. The cooling 
rate £ may be represented in the following way: 



c 



c* 
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(13) 



and Cjj may be obtained similarly [TB]. Here te is 
the Enskog relaxation time 



f(t) = -z92{v)(J 2 n 



nT(t) 



(14) 



2.3 Chapman-Enskog scheme 



To find f s we need to solve the Boltzmann-Lorentz 
equation (jfjj). It may be done approximately, using 
the Chapman-Enskog approach based on two simpli- 
fying assumptions: (i) f s (v, uj, r, t) depends on space 
and time only trough the macroscopic hydrodynamic 
fields and (ii) f s (v, uj, r, t) can be expanded in terms 
of the field gradients. In the gradient expansion 



= + Xf 



f(i) 



+ A 2 /.. 



(2) 



(15) 



where only leading terms in the expansion are kept. 



a formal parameter A is introduced, which indicates 
the power of the field gradient; at the end of the 
computations it is set to unity, A = 1. 

In a homogeneous cooling state with a lack of 
macroscopic fluxes only two hydrodynamic fields are 



3 



relevant - the number density of tracers, n s (r, t), and 
the total temperature T tot = T(t)+R(t). The former 
corresponds to conservation of particles, and the lat- 
ter to conservation of energy in the elastic limit. T tot 
is equal for the gas and tagged particles due to their 
mechanical identity. In the scaling regime, when the 
Chapman-Enskog approach is applicable, the trans- 
lational, rotational and total temperature are linearly 
related. Hence one can use any of these fields in the 
Chapman-Enskog scheme and we choose T(t) here 
without the loss of generality. Thus, the two rele- 
vant hydrodynamic fields satisfy in the homogeneous 
cooling state the following equations: 



dt 



— = I h A 



dt V dt 



d& 

dt +A ~m 



where d^ip/dt indicates that only terms of fc-th or- 
der with respect to the field gradients are accounted. 



Substituting Eq. (15 1 into the Boltzmann-Lorentz 
equation ^ and collecting terms of the same order 
of A (that is of the same order in gradients) we obtain 



successive equations for fg (v, C3,f,t), fs 1 ' (v, w, r , i) , 
etc. The zeroth-order equation in the gradients yields 
fs°\v,uj,r,t) for the homogeneous cooling state. 
Since the tracers are mechanically identical to the 
rest of the particles, fs°\v,uj,r,t) is simply propor- 
tional to the distribution function of the embedding 
gas: 



fi°Hv,u,r,t) = n ^ t lf(v,u},t). 
n 

The first-order equation then reads 



(17) 



fi Q) = 92(a)! (f,f^). (18) 



Using Eqs. (16), one can show, that the first term in 



the left-hand side of Eq. ( 18 ) is equal to 



dt 



dt dn. 



dt dT 



= -CT 



df, 



(i) 



dT 



(19) 



since d^°'n a /dt. Similarly we obtain that 
d^fP/dt = 0. Substituting Eq. (rr| with the space 
independent n and / into the last term in the left- 



hand side of Eq. ( 18 ), we arrive finally at the equation 
for/P: 



CT- 



df {1) / 



f 



v\ ■ Vii, 



(20) 



We search for the solution of Eq. ( 20 ) in the form 

/« =G(v 1 ,w 1 ,t)-Vn s (r,t) , (21) 



which implies with Eq. ( 10 1 the diffusion flux 



X 2 DV 2 n s 



T = -(T 



J s = I dvi / dw x vi [G ■ Vn s j = -~DVn s , (22) 
(16) 

where we take into account the isotropy o f ze ro-order 
function /i 0) = / s (0) (|t/|,|cj|). From Eq. Q we ob- 
tain (cf. [5]): 

(23) 



£) = ——/ ill-} I ilJ\ • Ct , 



We substitute / s (1) from Eq. |21j) into Eq. |20f dis- 
carding the factor Vn s . Then we multiply it by |ui 
and integrate over v\ and Q\ to obtain 



dv\ \ duj\ x 



Xv 1 -l(f,G)=^ m Jdv.jdcJ^f. (24) 

The first term in the left-hand side equals 
-(TdD/dT, while the right-hand side equals (T/m), 
according to Eq. ( lTj) . The structure of Eq. ( 20 1 dic- 



tates the Ansat£ 

3(v,fi,t)(xvf(v,L3,t)=B vf, (25) 

where we assume that the mean value of the rota- 
tional velocity to is equal to zero. The unknown con- 
stant B is to be determined from the above equation. 

2 Only the first term of the expansion B = B n + 
BlS^Cu 2 /■»§,) + B 2 sf^(v 2 /v*) + ... is used in Eq. see 
e.g. 
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With this Ansatz the second term in the left-hand 



side of Eq. ( 24 ) reads 

92 (<r) , ,_ ,_. 

aujidviVil 



(f,G) 



g 2 (<j)a 2 B 



where I?£;(t) = (T(i)/m) te(£) is the Enskog value 
for the self-diffusion coefficient for given temperature 
T(t). Note that for r\ — 0, i.e., for smooth particles, 



3 2 

x J diO\duj2 J de 9 (— U12 • e) \vi2 ■ e | (26) 
x / (vi,ui)f (V 2 ,w 2 ) (vi - v 2 ) ■ (v[ - Vi) . 



dv\dv2 Eq. (31 ) reproduces the previously known result, e.g. 



To evaluate the expression in Eq. ( 26 ) we use the col- 
lision rule ([5]) for the factor (v[ — v\) and Eq. ( f 2 ) for 
the distribution function /. After a straightforward 
algebra we arrive at 

32 [a) 



3 

where 



duji / dviV\I 



(f,G) 



. ad 



v 



1 



16 



B Tn _ x 

™ , ad ' 



a-20 



(27) 



(28) 



Note that t v a( j has a physical meaning of the relax- 
ation time of the velocity correlation function of gran- 



ular particles, e.g. 0(5]. Eq. (28 1 gives the general- 
ization of this quantity for the case of rough grains. 
If we express the unknown constant B$ using the re- 
lation 



dvi / dc3i«i • Vif(vi,Si,t) = - 



according to Eqs. (23 11), we recast (24 1 into the 
form 



,dD 
df 



T 



(29) 



As it follows from Eqs. (13), (14) and (28), T~{ d cx 
fiT and CT-j^p oc C oc cx -v/T- At the same 



oc 



time, the right-hand side of Eq. (29) scales as cx T, 
which implies that D cx y/T. Therefore T§f = f 



and we obtain the solution of Eq. ( 29 ) 
T 



£>(*) = 



1 



^Lw - «C(*) 



(30) 



Substituting r„ jac i and £ given by Eqs. (13) and (28) 



into Eq. (30) we finally find: 

D = 



Of 



(1+^20) 



(i- 



Q20 

16 



ail 

4 



(31) 



All the expansion coefficients 020, an, 002 and 6 
are small and with a reasonable accuracy one can use 
the Maxwell approximation for r [T2l [18] : 



= v 7 ! + C 2 + C 



(32) 



C 



2g(l + /3) 



l + (3 



(l + g)-(l-g)(l-j9) 



and for D, which simplifies in this case to 

-1 



D = D 



E 



+ 
4 



»7 



1 



(33) 



It is rather straightforward to perform similar cal- 
culations for a driven granular gas in the white-noise 
thermostat. The Boltzmann equation reads in this 
case: 



d x 2 



f = : 92 (a) /(/,/), (34) 



where \ characterizes the strength of the stochastic 
force. The gas temperature and the coefficients a 2 o, 
a iii a 22 are then determined by the intensity of the 
noise [18) . The temperature rapidly relaxes to a con- 
stant value, therefore only one hydrodynamic field 
n s (r,t) is relevant in the Chapman-Enskog scheme. 
Performing then all steps as in the case of HCS, we fi- 
nally arrive at the self-diffusion coefficient for a driven 
granular gas: 



Ah. 



' 1+g 

. 2 



(35) 



3 Molecular dynamics simula- 
tions 

To check the predictions of our theory we perform 
molecular dynamics (MD) simulations [T5] for a force- 
free system, using 8000 spherical particles of ra- 
dius a/2 = 1 and mass m = 1 in a box of length 
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£box = 132 with periodic boundary conditions. We 
confirmed that the system remained in the HCS dur- 
ing its evolution. The initial translational and rota- 
tional temperatures were equal T(0) = R(0) = 1; the 
ratio between translational and rotational tempera- 
tures r rapidly relaxed to a constant value, indicating 
that the application of our theory, based on the scal- 
ing form of the distribution function, is valid. To 
analyze the diffusion coefficient we used the re-scaled 
time r 



dr = dt/T c (t), r~ x {t) = A^g 2 {(j)(J 2 n^T{t)/m. 

(36) 

In the re-scaled time the self-diffusion coefficient at- 
tains at t > 1 a constant value [see Eq. Q], D = 
(R 2 (t))/(6t). 

The dependence of the self-diffusion coefficients on 
the normal and tangential restitution coefficients is 
presented in Fig. [I] where the molecular dynamics 
data are compared with the theoretical predictions. 
As it may be seen from the figure, the relative diffu- 
sion coefficient D(t)/DE(t) increases with decreasing 
e - the dependence, which has been already observed 
for smooth particles H El |2] . The physical na- 
ture of the effect is very simple: With the increasing 
inelasticity, which suppresses the normal component 
of the after-collisional relative velocity of particles, 
their trajectories become more stretched. This leads 
to the increasing correlation time and and hence, to 
a larger D. At the same time self-diffusion coefficient 
decreases for large roughness (/3 > —0.25) with in- 
creasing tangential restitution coefficient, see Fig. [T] 
It is not difficult to explain the observed behavior 
of D. Indeed, when the tangential restitution coeffi- 
cient j3 increases from —1 (smooth particles) to +1 
(absolutely rough particles), the translational and ro- 
tational motion become more and more engaged [see 
Eqs. ([5])] and the trajectories of particles more and 
more chaotic, that is, less stretched. This eventually 
causes the decrease of the diffusion coefficient with 
roughness. 

As it is follows from the Fig. [T] the theoretical 
predictions for the self-diffusion coefficient are in a 
reasonably good agreement with the numerical data. 
The diffusion coefficient, calculated in the framework 
of the Maxwellian approximation, practically does 
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Figure 1: The dependence of the self-diffusion coeffi- 
cient in a HCS on the normal e (upper panel) and tan- 
gential (3 (lower panel) restitution coefficients. The 
MD simulation data (symbols) are compared with the 
predictions of the theory for a gas in a HCS, Eq. (31 1 



(solid line) and Eq. (33) (dashed line, Maxwellian ap- 
proximation). The fixed restitution coefficients are: 
j3 = 0.9 for the upper and s = 0.9 for the lower panel. 
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not differ from the full solution, Eq. (33); the dif- References 



ference becomes apparent only for very high inelas- 
ticity (e < 0.6). A slight overestimate by the the- 
ory of the value of D, obtained in molecular dynam- 
ics, may be possibly attributed to the simple Ansatz 



(|25| for the first-order function fs , where only zero- 



order term of the expansion in the Sonine polyno- 
mials Sy\(v 2 /Vj) was used (see the footnote before 



Eq. (25)). One probably needs to extend the expan- 



sion and include the next order terms; this would be 
a subject of a future study. 
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